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Abstract 



Restoration of digital images from their degraded measurements has 
always been a problem of great theoretical and practical importance in nu- 
mcrous applications of imaging sciences. A specific solution to the problem 

U'^ of image restoration is generally determined by the nature of degradation 

phenomenon as well as by the statistical properties of measurement noises. 
C/J The present study is concerned with the case in which the images of intor- 

, est are corrupted by convolutional blurs and Poisson noises. To deal with 

such problems, there exists a range of solution methods which are based 
on the principles originating from the fixed-point algorithm of Richardson 
^ and Lucy (RL). In this paper, we provide conceptual and experimen- 

l/^ tal proof that such methods tend to converge to sparse solutions, which 

makes them applicable only to those images which can be represented 
by a relatively small number of non-zero samples in the spatial domain. 
Unfortunately, the set of such images is relatively small, which restricts 
the applicability of RL-type methods. On the other hand, virtually all 
practical images admit sparse representations in the domain of a properly 
CT^ designed linear transform. To take advantage of this fact, it is therefore 

tempting to modify the RL algorithm so as to make it recover represen- 
^ tation coefficients, rather than the values of their associated image. Such 

• ^ modification is introduced in this paper. Apart from the generality of its 

assumptions, the proposed method is also superior to many established 
^ reconstruction approaches in terms of estimation accuracy and computa- 

tional complexity. This and other conclusions of this study are validated 
through a series of numerical experiments. 



1 Introduction 

The notion of event counts is fundamental in many imaging modalities. Thus, 
for instance, this notion is used to quantify the reception of gamma photons in 
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nuclear imaging [l||4| as well as to describe the formation of optical images in 
charge coupled devices (CCD) [5j|6]. Confocal microscopy astronomical [sj 
and turbulent imaging [9] are additional examples of applications in which the 
notion of even counts routinely arises. In all these cases, working with event 
counts entails using a specific statistical assumption on the nature of acquired 
images. In particular, the images are assumed to be formed as blurred versions 
of their original counterparts contaminated by Poisson noise. Consequently, to 
recover the original images, the combined effect of the blur and noise needs to 
be reversed. A novel approach to the solution of this classical inverse problem 
is proposed in this work. 

Since in practice, digital images are discrete and finite-dimensional objects, 
it seems reasonable to represent the original image / along with its measured 
version g hy N x M real- valued matrices. Moreover, since in Poissonian imaging 
the values of both / and g have the interpretation of event counts, the images can 
be further constrained to be members of the subspace V of nonnegative-valued 
matrices. Formally, 

f,geY={ye M^^^ | y„.™ > 0, n = 0, 1, . . . , iV - 1, m = 0, 1, . . . , M - l} . 

(1) 

Additionally, let H : V — > V be a linear operator describing the convolution of 
/ with a non-negative mask ft, > of some size. Then, the formation of g can 
be formally expressed as 

g^vmm, (2) 

where V describes the effect of contamination of by Poisson noise. There- 

fore, as was mentioned before, to recover the original image / from g, the 
combined effect of T'j'Hl-}} in ([2| needs to be inverted. 

The current arsenal of image restoration approaches to the solution of ([2| 
is relatively broad. In general, all these methods can be divided into two main 
groups, where the methods of the first group assume Ji to be identity, whereas 
the methods of the second group allow it to be a general low-pass filter. Within 
the first group, a range of available reconstruction methods take advantage 
of a variance stabilized transformation (VST) [10| , which allows the Poisson 
noise in ^ to be transformed into approximately white Gaussian noise. Some 



recent developments in this direction include the works in 11 ■ 13 , in which 



wavelet de-noising has been exploited to reject the "gaussianized" Poisson noise. 



Conceptually similar ideas have been also advocated in 14 15 based on more 



advanced tools of statistical analysis. The framework of Bayesian estimation was 



exploited in 16 -20 , where hierarchical (Markovian) models have been used to 
describe a priori probabilities of the original images. Although all the above 
methods can be used to recover an approximation of / in ([2]) , their applicability 
is restricted to the case of weak blurs, where can be closely approximated by 
the identity operator Id. Moreover, all these methods depend on VST, whose 



performance is known to deteriorate considerably in low-count settings 21 . 

One of the nowadays classical methods of solving ([2| for the case of a gen- 
eral H was proposed in the beginning of the 1970s in the seminal works of 
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W. Richardson 22 and L. Luci 23 . This method - known as the Richardson- 
Lucy (RL) algorithm - can be classified as a maximum likelihood (ML) esti- 
mator. As a general rule, however, ML estimation may result in less accurate 
and/or stable reconstructions as compared with maximum a posteriori (MAP) 
estimation methods based on the Bayesian paradigm. Consequently, the RL 
algorithm has been recently extended under the MAP estimation framework, 
resulting in a number of regularized solutions to ^ which differ in the way the 
original image / is modelled. Thus, for example, Gaussian models have been 
exploited in 24 25 , while the algorithm in [26] is based on the total variation 



model of 27 . Unfortunately, there are conditions on which the above methods 
can result in erroneous reconstructions (as will be demonstrated later in the 
paper) . 

Despite the relative simplicity of the RL algorithm, it still remains one of 
the most widespread methods used in the current practice of Poissonian imag- 
ing. The main advantage of the RL algorithm stems from its nonlinear nature, 
allowing one to recover the high-frequency components of the original images 
which are the most affected by blur. Moreover, a closer look at the analytical 
properties of the RL algorithm reveals the fact that it has a "built-in" ability to 
converge to sparse solutions, and, as a result, the application of the algorithm 
for the restoration of sparse images should be expected to produce particularly 
useful reconstructions. Unfortunately, most of the real-life images are not sparse 
in the spatial domain, in which case the RL algorithm may fail to provide useful 
results. On the other hand, most of such images are sparse in the domain of a 



properly designed linear transform 28 -31 . Consequently, to extend the appli 



cability of the RL algorithm to a wider range of imagery data, it is tempting to 
find a way to apply the algorithm in the transformed domain, as opposed to the 
spatial domain. Accordingly, the main contribution of this paper consists in the 
introduction of a novel approach to the solution of ([2| which exploits the above 
idea. Moreover, it will be shown via an experimental study that the proposed 
method outperforms a number of alternative algorithms in terms of normalized 



mean-square error (NMSE), SSIM quality index 32 , as well as its stability and 
computational efficiency. 

The rest of the paper is organized as follows. Section II provides some 
necessary details on the RL algorithm. The statistical models and assumptions 
underpinning the proposed method are detailed in Section III, while Section IV 
summarizes the main structure of the proposed solution method. The results of 
both computed-simulated and real-life experiments are presented in Section V. 
Section VI finalizes the paper with a discussion and conclusions. 



2 Richardson-Lucy Approach 

To establish a necessary foundation for the development and presentation of the 
proposed method, a brief overview of the RL algorithm is provided first. Under 
the image formation model of (l2| , the likelihood function of the measured image 
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g can be defined as 



Af-l M-1 -(«{/})„ „ f-Tjr f 1 xS„,™ 

^(g I /) - n n ^^^^""'"^ - (3) 

71—0 m— ' 

where the samples ^n.m of 9 are assumed to be mutuahy independent. Con- 
sequently, defining E{f) = — logP(^ | /), the ML estimate /ml of / is given 

by 

Iml = argmin{£;(/)}, (4) 

/ev 

where 

E{f) = {l,n{f})-{g,logmf})), (5) 

with (•, •) denoting the standard inner product in Differentiating E{f) 

with respect to / and equating the derivative to zero results in 

where 1 denotes an x M matrix of ones, the fractional line stands for a 
point-wise division, and H* : V — > V denotes the adjoint of H. Note that if H 
represents the operation of 2-D convolution with a positive- valued kernel h > 0, 
then H* represents the convolution with a spatially reversed version of the same 
kernel h. Moreover, if h is normalized to satisfy ^„ ^ hn,m = 1, then both 
and T-1*{1} are obviously equal to 1. This fact is exploited by the RL algorithm, 
which recovers an approximation of the original image / as a stationary point 
of the sequence of solutions produced by 

f (t+i) ^ fit) . ^* / 9 1 z^-, 

where the dot stands for a point-wise multiplication, and t is the iteration index. 

To better understand the nature of the solution produced by the minimiza- 
tion of ([5]), consider the following steps. Using the definition of an adjoint 
operator, the cost functional E{f) in ^ can be rewritten as 

E{f) = {n*{l}J)-{g,\og{n{m= (8) 
= (l,/)-(5,log(H{/})), 

Moreover, since / is assumed to be positive valued, the inner product (1, /) 
coincides with the ^i-norm of / defined as ||/||i = J2n m \ fn,m\- As a result, the 
cost function E{f) in Q can be further redefined as 

£^(/)-||/lli-(5,log(H{/})). (9) 

A closer look into the two terms in (|9| leads to a number of remarkable 
observations. First and foremost, the presence of ^i-norm in the cost function 
implies that the solution of Q will tend to be sparse 28 29 33 - 35 . Therefo 



fore, 
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the RL algorithm has a "built-in" ability to converge to sparse solutions. More- 
over, such solutions are guaranteed to be positive due to the second term in 
([9| which works akin to a logarithmic barrier [36| Ch.ll]. This is what makes 
the RL algorithm to favour the reconstructions which are sparse in the spatial 
domain. 

It goes without saying that the assumption of spatially sparse images may be 
unacceptable for a variety of natural scenes. For this reason, the RL algorithm is 



known to perform reliably in the case of, e.g., astronomical images 22 , while its 



application to piecewise smooth images can produce rather disappointing results 



26 . On the other hand, sparsity is known to be an extremely useful and liable 
assumption to use for describing the behaviour of the representation coefficients 
of natural images in the domain of certain linear transforms. Therefore, to 
extend the applicability of the RL algorithm to the above case, it is imperative 
to find a way to apply the algorithm directly to the (sparse) representation 
coefficients, rather than to their corresponding images. A possible variant of 
such an approach is detailed in the section that follows. 

Before turning to the description of our method, a number of important 
properties of the RL recursion in Q are worth to be mentioned. In particular, 
it can be seen from ([t]) that the algorithm preserves the positivity of /^*\ viz. 
j(t+i) jg guaranteed to be in V provided Z*^*^ G V and g € Y. Somewhat less 
trivial, it can also be shown that the algorithm preserves the mean values of the 
intermediate solutions, i.e. ^„ ,„5n,m = J2n m fnln at each iteration t 7 . All 
these properties will play an important role in the proposed method as detailed 
below. 



3 MAP Formulation 

The approach proposed in this paper is based on the assumption that the original 
image / in ([2| can be sparsely represented in the domain of a linear transforma- 
tion. In particular, given an overcomplete and dense set {ipklkei of vectors in 
jgjVxM^ we define the synthesis operator $ to be given by 

*:£2(2:)^V:cK^*{c} = ^Cfc0fe, (10) 

fcei 

where T denotes the set of indices of the representation coefficients c, with 
NM < jfX < oo. Using this definition of one can rewrite the model equation 
([2| as 

.9 nA{c}}, (11) 

with A being a composition of H and i.e. A := 7i o ^, and c e ^2(2^) being 
the representation coefficients of / in the domain of In this new format, the 
original problem of reconstruction of / is replaced by the problem of estimating 
its representation coefficients c. 

Needless to say, due to the over-completeness of the definition of c is 
not unique. This difficultly, however, can be easily overcome by requiring the 
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representation coefficients c to be as sparse as possible. A possible way to quan- 
tify the sparseness of c is by measuring its ii norm ||c||i 29 . It is worthwhile 



noting that this measure of sparseness has been successfully used in numerous 



fields of science, which include signal modeling 29 , compressed sensing [34[|35| , 
independent component analysis and blind source separation 37 38 , inverse 



39 


40 


sis 
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, signal and image de-noising 28 30 , morphological compo- 
together with its earlier version introduced in 42 . 
Using the assumptions above, the problem of recovering the original image 
/ can be cast into the framework of MAP estimation, in which case the optimal 
/map is found as a solution of the following maximization problem 



Imap = argmax{P(5r | c) P(c)} , 



(12) 



where P{g \ c) and P{c) denote the likelihood function and the prior probability 
of c, respectively. 

Under the assumption of statistical independence, the likelihood P{g \ c) 
has the form of 



Pig I c) 



nn — 



n.m 



n—0 m— 



(13) 



One the other hand, congruent to the assumption on minimality of |lc|li, the 
representation coefficients are assumed to be i.i.d. random variables that obey 
a zero-mean Laplacian distribution, viz 



^(-^^ = n 



(14) 



where 7 > is a scale parameter of the distribution. 



Using the definitions ( 13 1 and ( 14 1 and applying the standard log-transformation 



and the negation of sign to (12 1, it is straightforward to show that the MAP 



estimation problem is equivalent to 

Imap = argmin{£;(c)}, 

C 

Eic)^{l,A{c})-{g,logiA{c})) 



M\c\\ 



(15) 
(16) 



with A = 1/7 being a regularization parameter. 

Care should be exercised when specifying the domain of definition of E{c) in 
( 16 ). Indeed, the likelihood model of ( 13 ) interprets the value fn,m as the mean 
of the corresponding random observation gn,m- Since, in the case of a Poisson 
distribution, its first and second moments are equal, the values fn,m should be 
assumed to be non-negative (in accordance with our earlier assumption in ([T 
Thus, the domain of the objective E{c) is formally given by 



dom£; = {c 6^2(2:) I (*{c})„,™ > 0,Vn,m}, 



(17) 



6 



which is a convex set. Moreover, as long as gn,m > and the convolution 
operator H is non-degenerate (albeit, possibly, ill-conditioned), the objective 
functional E{c) is guaranteed to be strictly convex. In this case, E{c) admits 
a unique minimizer in domi?, which can be found by any optimization al- 
gorithm which is guaranteed to converge to a stationary point of E{c) [36] . 
Unfortunately, when either g is not strictly positive or H possesses a non-trivial 
null-space, the convexity of E{c) is not strict, and, as a result, the existence and 
uniqueness of its global minimizer cannot be a priori guaranteed. 



4 Solution via Fixed-Point Iterations 

Numerous algorithms approaches can be used to minimize the cost functional 



in (16 1. In this paper, we propose a different method based on fixed-point 



iterations. To this end, we first notice that the first-order optimality condition 



for ( 16 1 has the following form 



VE{c - c*) A*{1} - A* [jf^j + Asign(c*) = 0, (18) 
where A* denotes the adjoint of A and c* is an extremum of E{c). Note that, in 



(18), the non-differentiability of the absolute value at zero is resolved by letting 
(|x|)'|^_g — 0. As will be shown below, the latter assumption is rather technical, 
as it has no impact on the proposed solution. 



Finally, by rearranging the terms in (18), one obtains 



^*|^|=^*{1} + Asign(c*), (19) 

which, in turn, suggests the following fixed-point iteration algorithm 

r(*+i) - A* I - 1 — (20) 

M{c(*)}/ ^*{1} + Asign(c(*))- ^ ' 



The iterative procedure of ( 20 ) is multiplicative in nature - the fact that has 



a number of important implications. First, the zero entries of c are preserved 
by the iteration procedure, which means that if c|*^ = for some i e I, then 

c\ — for all t' > t. In this respect, letting (|a;|)' — sign(a;) seems to be a 
reasonable simplification, since the value of the gradient of ||c||i at zero does not 



seem to have any influence on the result of the iterative procedure (20). Second, 
if the values of c were allowed to be of an arbitrary sign, it would be generally 
impossible to guarantee a monotone convergence ofc(*). To overcome this de- 
ficiency, in what follows the representation coefficients c will be assumed to be 
non-negative. Moreover, if the representation vectors {(/>i}igx are constrained 
to be positive-valued as well, it is straightforward to verify that the iterations 



in ( 20 ) will preserve the positiveness of c. It is important to note that the above 



constraints on the sign of c and {0i}igx should not be seen as a limitation, since 
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there exists a body of works (see, for example, [43]), which describe the con- 
struction of positive- valued representation sets that allows the natural scenes to 
be represented in terms of both sparse and positive coefficients. 

Provided the elements of {4>i]i£i are positive values, the vector v :— A*{1\ 



{l]n,m('/'i)n,m| wiU be positivc valucd as weU, i.e. Vi 
can be rewritten in a more compact form as 



> 0,Vi. In this case, (201 



„(t+i) 



A* 



4t) 



(21) 



The iterations in (21) can be initialized with a constant coefficient vector c^^^ 



(e.g. c- = l,Vi) and executed until the relative change ||c^*+^ — c*^*' ||2/ ||c^'' ||2 
drops below a predefined threshold < e ^ 1. The above algorithm will be 
referred below to as the sparse RL method, or simply SRL. 

Comparing the recursions in (21 1 and ([7|, one cannot avoid noting how 
similar the SRL and RL algorithms are. Indeed, replacing c and A in (l2T|) by 



/ and respectively, yields an update equation identical to that in ([7| (up to 
the element-wise division by w + A in (21)). Moreover, the assumption on the 
non-negativeness of c and A made by SRL is parallel to the non-negativeness 
of / and T-L exploited by the RL method. Further similarity between RL and 
SRL reveals itself in the objective functions of the two reconstruction methods. 



Specifically, the objective function in ( 16 ) can be rewritten as 



E{c) - {A*{l},c) - (g,log(^{c})) + A||ci|i = 
= {v,c) - (g,log(^{c})) +A||c||i = ||c||^,: 



(.g,log(^{c})). 



(22) 



where w = v + X and ||c||tu,i = X^iei^i'^i w-weighted ^i-norm of the 

non-negative representation coefficients c. Thus, one can see that under the 
substitution of c by /, by "H and w — \, the objective functionals (22 ) and (|9| 
are identical. Yet, while minimizing ([9| forces the reconstruction of / to have 
a sparse appearance, minimizing (22) has the same effect on the representation 



coefficients. As opposed to the former, the latter case allows one to recover much 
broader classes of practical images, as it is demonstrated by the experimental 
results below. 



5 Experimental Results 

5.1 Reference methods and comparison measures 

Both one- and two-dimensional data sets have been used to test the performance 
of the proposed SRL algorithm. In the case of 1-D data, the SRL algorithm was 
compared against the RL method, while in the case of 2-D data, the comparison 



was made against the method described in 26 . Note that, just like SRL, the 
reference method of |26j belongs to the class of Bayesian estimators. Specifically, 
the method recovers the original image / as a maximizer of its posterior proba- 
bility, computed under the model of ^ and the total-variation (TV) prior. The 
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maximization can be performed itcratively according to the following update 
equation 

^*+'=l-7div(vVl|V/.||) "^infe}}' ^'^^ 

where 7 > is a regularization parameter, which has been set to be equal 
to 0.002 as prescribed in [2^. For the convenience of referencing, the above 
algorithm will be referred below to as the RLTV method. 

To quantitatively compare between the performance of various reconstruc- 
tion algorithms, two comparison measures have been used. The first, normalized 
mean-squared error (NMSE), is defined as follows. Let / be an original image 
and / be an estimate of /. Then, the NMSE can be defined as 



11/ -/I 



2 



NMSE = £ { \ , (24) 



with II • \\f being the Frobenius matrix norm, and £ being the operator of 
expectation. In the current study, the expectation was approximated by its 
corresponding sample mean based on the results of 200 independent trials. 

It has been recently argued that the NMSE may not be an optimal compar- 
ison measure as long as human visual perception is concerned. For this reason. 



the structural similarity (SSIM) index proposed by 32 was employed as well as 
an alternative performance metric. 



5.2 One-dimensional Reconstruction 

The first part of our experimental study is concerned with the recovery of piece- 
wise constant signals of length N. For this case, a natural choice would be to 
define the basis functions to be scaled and spatially shifted versions 

of a rectangular (Haar) box. Specifically, let j = 0, 1, ... , [log2(iV)J — 1 be a 
(dyadic) resolution index, and (f)j e be defined as 




for < n < 2^ - 1 
for 2J <n<N, 



(25) 



where the normalization factor 2 ■'/^ is used to guarantee ||0j||2 = 1 for all j. 
Also, let 

Z''^ : -> : x[n] ^ x[{n - k) modiV] (26) 



be the operator of (causal) circular shift by k points. Obviously, for each (pj there 
are a total of N non-repetitive shifts. Consistent with the notations standard 
in wavelet theory, let 0^ ^ := Z^{(j)j} denote the function (f)j (circularly) shifted 
by k points to the right (in which case (pk^j can be viewed as a rectangular box 
function supported on [{k, k + 2^ — l)modA^]). Moreover, let ^ := {(j>k,j} be 
a matrix of size N x iV[log2(A^)J , whose columns are formed by functions (j>k,j 
with A; = 0, 1, . . . , TV — 1 and j = 0, 1, . . . , J — 1. In this case, given an arbitrary 



9 



vector of representation coefficients c of size A'^[log2(A^)J x 1, its corresponding 
signal / can be synthesized as / = $ c. 

The dictionary $ constructed above is severely overcomplete, which has the 
disadvantage of high computational complexity associated with the computa- 
tions of # c and $^/. To improve the computational efficiency, the index j was 
restricted to four resolution levels, viz. j — 2, 3, 4, and 5. As a result, the size of 
$ was equal to 128 x 512 for N = 128. The coefficients c used for the synthesis 
of test signals had about 1.5 3 % of non-zero entries, drawn from a uniform 
distribution. The blurring artifact was simulated by convolving the test images 
with a Gaussian kernel H whose -3 dB cut-off frequency was set to be equal 
to 0.27r. The maxima of the resulting signals were set to two different values, 
namely 256 and 32, in order to simulate high- and low-count detection scenarios, 
respectively. As a final step, the blurred signals were contaminated by Poisson 
noise (see the second subplots in Fig. [T] and Fig. [2] for typical examples of data 
signals) . 

Next, the RL and SRL algorithms were applied to the synthesized data 
signals with A = 0.2. Unfortunately, no monotonous convergence in NMSE 
could be achieved in the case of the RL algorithm. The steady-sate estimates 
obtained using this classical method were unacceptably noisy (see the fourth 
subplots of Fig. [l] and Fig. [2] for typical examples of the RL estimates) . For 
this reason, it was decided to terminate the RL estimations before the steady- 
state was reached, at the point when the NMSE value was minimized. (Note 
that such NMSE-optimal RL estimates are impossible to compute in real-life 
scenarios, as their computation requires knowing the original images.) The SRL 
algorithm, on the other hand, produced monotonous convergence in NMSE, 
producing substantially smaller values of NMSE as compared to the case of 
RL estimation. For the case of high- and low-count scenarios, some typical 
reconstruction results are depicted in Fig. [l] and Fig. [2] respectively, from which 
the superiority of the proposed method can be easily appreciated. It is also 
interesting to note that, as predicted by the theory, the RL algorithm attempts 
to arrive at a sparse estimate of the signal of interest, which contradicts the 
piece-wise smooth nature of our test signals. The SRL method, on the other 
hand, require the sparsity of the representation coefficients, which appears to 
be a much more realistic assumption for the case at hand. 

A quantitative comparison of the two reconstruction algorithms is presented 
in Fig.[3j which depicts the values of NMSE produced by the RL and SRL meth- 
ods as a function of the number of iterations. It is important to emphasize that 
each value of the NMSE in Fig. [3] is a result of averaging the errors obtained 
in a series of independent trials, where both the true signals and noises were 
drawn randomly. Observing Fig. [3] one can see that SRL results in considerably 
lower values of the NMSE as compared to the case of RL. As well, it converges 
monotonously to a steady-state solution after a relatively small number of it- 
erations. On the other hand, the convergence of RL is not monotone ~ the 
fact that necessitates termination of the algorithm after a predefined number 
of iterations. Needless to say, determining a proper number of iterations is a 
data-dependent and somewhat esoteric task in practice. 
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MEASURED SIGNAL, NMSE = 0.03 




RL ESTIMATE (NMSE-OPTIMAL), NMSE = 0.027 




RL ESTIMATE (STEADY STATE), NMSE = 0.433 



SPARSE RL ESTIMATE, NMSE = 0.002 




Figure 1: High-count scenario (from up to bottom): The original signal, its 

blurred version, the blurred signal contaminated by Poisson noise, the MSE 
optimal RL estimate, the steady-state RL estimate, and the SRL estimate. 




Figure 2: Low-count scenario (from up to bottom): The original signal, its 
blurred version, the blurred signal contaminated by Poisson noise, the MSE 
optimal RL estimate, the steady-state RL estimate, and the SRL estimate. 
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Figure 3: The convergence of NMSE values produced by the RL and SRL 
algorithms in the case of high (upper subplot) and low (lower subplot) count 
scenarios. 
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5.3 Two-dimensional Reconstruction 



In the second part of the experimental study, the SRL algorithm is examined 
using imaging data. As reference methods, the standard RL method as well as 
the RLTV algorithm (as specified by (23)) are employed to assess the perfor- 
mance of SRL in a comparative way. Note that as neither RL nor RLTV could 
provide acceptable steady-state results, their iterations needed to be terminated 
at an earlier point. In particular, the iterations were terminated at the point 
when the NMSE produced by the algorithms reached a minimum value. It is 
important to reiterate that such "MSE-optimal" termination rule is not realiz- 
able in real-life scenarios where the original images are unknown. Hence, the 
RL and RLTV estimates demonstrated below represent the best possible result 
which could be achieved using these reconstruction methods. 

Implementation of the SRL algorithm relies on the availability of a dictio- 
nary of positive-valued atoms {4>i\iei- the 2-D experiments, we explore two 
possible approaches to the definition of such a dictionary, viz. unsupervised and 
supervised. In the first case, the dictionary is composed of shift-invariant sub- 



sets generated by cubic B-splines 44 , and therefore it does not depend on the 
nature of imagery data at hand. In the second case, the dictionary is learned 
from a training set that is supposed to represent the family of objects which the 
measured image is likely to belong to. 



5.3.1 Unsupervised dictionary 

For the sake of reproducibility of the results of this paper, we provide some 
necessary details on the construction of the cubic spline dictionary. To this end, 
let 12^x1 denote a vector of length 2-' , whose elements are all equal to 1. Also, 
let 63 = [1, 4, l]/6 be the vector of integral values of the cubic B-spline and hj 
be the vector of integral values of the cubic B-spline scaled by a dyadic factor 
of 2^ with j = 0, 1, . . . , J - 1. The vectors bj have the length of Mj- = 2^+^ - 1 
points and they can be computed recursively according to [44, 

bj = 2-3i ( l2^^^^.^^l2^ ) * 63 , (27) 

4 times 

where * denotes the operation of convolution. As prescribed by common prac- 
tice, the vectors hj can be subsequently normalized to obey \\bj\\2 = 1- 

The 1-D kernels bj can be extended into isotropic and separable 2-D (dis- 
crete) splines Bj by means of the tensor product, i.e. Bj[n,m] = bj[n]hj[m], 
< j < J —I. Consequently, the true image / is assumed to belong to the space 
of all integer translations of the discrete kernels Bj . Under periodic boundary 
conditions, there are a total of N x M representation coefficients Cj g E:^J^^*^ 
for each j = 0,l,...,J— 1. Let the array of all these coefficients be denoted by 
= {cjl/^o • Then, the operator A can be shown to be given by 

.7-1 

^ . ^WxMxJ ^ ^NxM . c J ^ ^1 ^ @ B,}, (28) 

3=0 
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while its adjoint A* is given by 
A 



!,NxM 



NxMxJ 



f^c={H{f}®B,} 



j-i 



(29) 



with ® standing for circular convolutiorQ According to ( 29 ) , an x M image 
produces a total oi J N M spline coefficients, and hence the overall complexity 
of applying A and A* is comparable to that of a stationary wavelet transform. 



5.3.2 Supervised dictionary 

In the supervised case, the dictionary {4>i}i£x is designed adaptively based on 
a set of training examples that represent the family of images, to which / is 
expected to belong. Consequently, the training procedure aims to find a set 
of positive-valued representation functions, using which one can represent the 
training images with a minimum possible number of representation coefficients. 
In this work, the training was performed by means of the the non-negative kSVD 
(NN-kSVD) algorithnj^ proposed in 43 . In compliance with the Sparseland 



model of [30], the training was applied to 16 x 16 overlapping segments of a 
total of 11 MRI scans of the brain, which did not include the scans used in 
reconstruction. Each of these segments was represented by a linear combination 
of a total of 512 atoms, thereby resulting in the over-completeness factor of 2 
per segment. For the formal definition of operators A and A* the reader is 



kindly referreded to 30 . 



5.3.3 Reconstruction of MRI scans 

The blurring kernel used in our next experiment was defined to be = 
-|- + 1)""^, with i,j — —7, . . . , 7 [31], while the regularization parameter 
A was set to be equal to 0.1 in both unsupervised and supervised cases. A 
typical example of the original image used in this study along with its blurred 
and noise-contaminated versions are shown in the leftmost, middle, and the 
rightmost subplots of Fig. |4] respectively. The "measured" MRI scans have 
been obtained by contaminating the blurred image with Poisson noise giving 
rise to an SNR of approximately 15 dB. 

Some typical reconstruction results produced by the proposed and reference 
methods are shown in Fig. [5j One can see that the SRL method is capable of 
better recovering the details of the true image as compared to the alternative 
approaches. It can also be seen that the SRL estimates possess higher resolution 
and contrast as compared to the alternative approaches. These conclusions are 
further supported by the quantitative measures of Table [T] which compares the 
NMSE and SSIM indices of different estimates. As shown by the table, the 
supervised SRL method produces the lowest NMSE and the largest SSIM index 
among all the methods under comparison. 

^For some examples of related MATLAB® codes and their use visit 
www. ece .uwaterloo . ca/'olegm/research.html. 

^A software for the NN-kSVD algorithm can be obtained from 
Ihttp : //www . cs . technion . ac . il/Telad/sof tware/ . 
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ORIGINAL IMAGE BLURRED IMAGE BLURRED & NOISY IMAGE 




Figure 4: (Left subplot) Original image of a cross section of the brain; (Middle 
subplot) Blurred image; (Right subplot) Blurred and noisy image. 



Table 1: NMSE and SSIM values obtained using the reconstruction methods 
under comparison. 





RL (MSE-optimal) 


RLTV (MSE-optimal) 


SRL (SpUnes) 


SRL (NN-kSVD) 


NMSE 


0.022 


0.011 


0.007 


0.003 


SIM 


0.71 


0.77 


0.89 


0.91 



6 Discussion and Conclusions 

In the current paper, a new regularized version of the RL algorithm has been pre- 
sented, which is advantageous in a mimbcr of ways. First, the proposed method 
is general in its formulation. The latter allows applying the same reconstruc- 
tion procedure to a number of different settings, such as image de-noising or 
image enhancement through dcconvohition. Attuned to deal with Poissonian 
noises, the method can therefore be applied for the reconstruction of imagery 
data produced by many important image modalities, including optical, micro- 
scopic, turbulent, and nuclear imaging, just to name a few. Moreover, whilst 
many alternative reconstruction methods take advantage of certain simplifying 
assumptions about the noise nature, the proposed technique is optimized to deal 
with the realistic noise model at hand. 

Another important advantage of the proposed method consists in the gener- 
ality of the prior model used to describe the nature of true images. Specifically, 
the images have been assumed to admit a sparse representation in the domain 
of a properly chosen linear transform. Note that the above a priori model is 
nowadays considered to be a standard model in numerous applications of the 
theory of sparse representations. This is what allows the SRL algorithm to re- 
cover piecewise smooth images - the task impossible to accomplish by means 
of the standard RL procedure for its property to favour spatially sparse recon- 
structions. 
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RL (MSE optimal) RLTV (MSE-optimal) 




SRL (NN-kSVD) SRL (splines) 




Figure 5: Image reconstruction results corresponding to Fig[4j (Upper row of 
subplots) The MSE-optimal RL and RLTV estimates; (Lower row of subplots) 
supervised (kSVD) and imsupervised (splines) SRL recovery. 
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Yet another critical advantage of the proposed SRL algorithm is in its al- 
gorithmic structure, which allows solving non-smooth optimization problems 
at the computational cost of a steepest descent procedure. Consequently, the 
computational load required by the proposed method is relatively small, which 
allows the method to be applied for the solution of large-scale problems and/or 
for processing of large data sets. 
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